% By Martin M. Andreasen

function [model,errorMes] = perturbationDSGEmodel(params,orderApp,setupModel,MexOn,selectY)

% Computing the steady state and checking that nf = 0
[auxOut,nf,errorMes] = numF(params);
params = auxOut.params;
if any(abs(nf)> 1D-8) 
    disp('Displaying nf:')
    disp(nf)
end
if errorMes ~= 0
    model = NaN;
    errorMes = 1;
    return;
end

% The dimension of the model
eta = auxOut.eta;
ne  = auxOut.ne;
nx  = auxOut.nx;
ny  = auxOut.ny;

% Ensure same ordering of params as in symparams
paramsOrdered = struct();
for i=1:length(setupModel.symparams)
    paramsOrdered.(char(setupModel.symparams(1,i))) = params.(char(setupModel.symparams(1,i)));
end
     
% Choose solution algorithm: 
% algo='vectorize'; % standard vectorization - good only for small models
%algo='dlyap'; % Hessenberg-Schur algorithm - good for large models, but
%                 requires to have the control system toolbox
algo='gensylv'; % Kamenik (2005) algorithm - recommended. If the MEX file gensylv.mexw64 does not work on your system, download dynare from www.dynare.org, search for gensylv and add the MEX file to the Perturbation folder.

% Solve the perturbation loadings
derivs = solve_dsge(setupModel,struct2array(paramsOrdered),...
    auxOut.momEps,eta,auxOut.xssTrans',auxOut.yssTrans',orderApp,algo,MexOn);


%% Saving the output in model
model.labelx           = auxOut.labelx;
model.labely           = auxOut.labely;
model.g0               = auxOut.yssTrans;
model.h0               = auxOut.xssTrans;

model.derivs           = derivs;
model.gx               = derivs.gx(1:ny,1:nx);
model.hx               = derivs.hx(1:nx,1:nx);
if orderApp > 1
    tmp_gxx            = reshape(derivs.gxx,ny,nx+1,nx+1);
    model.gxx          = tmp_gxx(1:ny,1:nx,1:nx);
    model.gss          = tmp_gxx(1:ny,nx+1,nx+1);
    tmp_hxx            = reshape(derivs.hxx,nx,nx+1,nx+1);
    model.hxx          = tmp_hxx(1:nx,1:nx,1:nx);
    model.hss          = tmp_hxx(1:nx,nx+1,nx+1);
else
    model.gxx          = zeros(ny,nx,nx);
    model.gss          = zeros(ny,1);
    model.hxx          = zeros(nx,nx,nx);
    model.hss          = zeros(nx,1);
end

if orderApp > 2
    tmp_gxxx           = reshape(derivs.gxxx,ny,nx+1,nx+1,nx+1);
    model.gxxx         = tmp_gxxx(1:ny,1:nx,1:nx,1:nx);
    model.gssx         = squeeze(tmp_gxxx(1:ny,nx+1,nx+1,1:nx));
    model.gsss         = squeeze(tmp_gxxx(1:ny,nx+1,nx+1,nx+1));
    tmp_hxxx           = reshape(derivs.hxxx,nx,nx+1,nx+1,nx+1);
    model.hxxx         = tmp_hxxx(1:nx,1:nx,1:nx,1:nx);
    model.hssx         = squeeze(tmp_hxxx(1:nx,nx+1,nx+1,1:nx));
    model.hsss         = squeeze(tmp_hxxx(1:nx,nx+1,nx+1,nx+1));
else
    model.gxxx         = zeros(ny,nx,nx,nx);
    model.gssx         = zeros(ny,nx);
    model.hxxx         = zeros(nx,nx,nx,nx);
    model.hssx         = zeros(nx,nx);
end
    
if orderApp > 3
    tmp_g4x            = reshape(derivs.gxxxx,ny,nx+1,nx+1,nx+1,nx+1);
    model.g4x          = tmp_g4x(1:ny,1:nx,1:nx,1:nx,1:nx);
    model.gssxx        = squeeze(tmp_g4x(1:ny,nx+1,nx+1,1:nx,1:nx));
    model.gsssx        = squeeze(tmp_g4x(1:ny,nx+1,nx+1,nx+1,1:nx));
    model.g4s          = squeeze(tmp_g4x(1:ny,nx+1,nx+1,nx+1,nx+1));
    tmp_h4x            = reshape(derivs.hxxxx,nx,nx+1,nx+1,nx+1,nx+1);
    model.h4x          = tmp_h4x(1:nx,1:nx,1:nx,1:nx,1:nx);
    model.hssxx        = squeeze(tmp_h4x(1:nx,nx+1,nx+1,1:nx,1:nx));
    model.hsssx        = squeeze(tmp_h4x(1:nx,nx+1,nx+1,nx+1,1:nx));
    model.h4s          = squeeze(tmp_h4x(1:nx,nx+1,nx+1,nx+1,nx+1));
else
    model.g4x          = zeros(ny,nx,nx,nx,nx);
    model.gssxx        = zeros(ny,nx,nx);
    model.gsssx        = zeros(ny,nx);
    model.g4s          = zeros(ny,1);
    model.h4x          = zeros(nx,nx,nx,nx,nx);
    model.hssxx        = zeros(nx,nx,nx);
    model.hsssx        = zeros(nx,nx);
    model.h4s          = zeros(nx,1);
end   
if orderApp > 4
    tmp_g5x            = reshape(derivs.gxxxxx,ny,nx+1,nx+1,nx+1,nx+1,nx+1);
    model.g5x          = squeeze(tmp_g5x(1:ny,1:nx,1:nx,1:nx,1:nx,1:nx));
    model.gssxxx       = squeeze(tmp_g5x(1:ny,nx+1,nx+1,1:nx,1:nx,1:nx));
    model.gsssxx       = squeeze(tmp_g5x(1:ny,nx+1,nx+1,nx+1,1:nx,1:nx));
    model.g4sx         = squeeze(tmp_g5x(1:ny,nx+1,nx+1,nx+1,nx+1,1:nx));
    model.g5s          = squeeze(tmp_g5x(1:ny,nx+1,nx+1,nx+1,nx+1,nx+1));
    tmp_h5x            = reshape(derivs.hxxxxx,nx,nx+1,nx+1,nx+1,nx+1,nx+1);
    model.h5x          = squeeze(tmp_h5x(1:nx,1:nx,1:nx,1:nx,1:nx,1:nx));
    model.hssxxx       = squeeze(tmp_h5x(1:nx,nx+1,nx+1,1:nx,1:nx,1:nx));
    model.hsssxx       = squeeze(tmp_h5x(1:nx,nx+1,nx+1,nx+1,1:nx,1:nx));
    model.h4sx         = squeeze(tmp_h5x(1:nx,nx+1,nx+1,nx+1,nx+1,1:nx));
    model.h5s          = squeeze(tmp_h5x(1:nx,nx+1,nx+1,nx+1,nx+1,nx+1));
else
    model.g5x          = zeros(ny,nx,nx,nx,nx,nx);
    model.gssxxx       = zeros(ny,nx,nx,nx);
    model.gsssxx       = zeros(ny,nx,nx);
    model.g4sx         = zeros(ny,nx);
    model.g5s          = zeros(ny,1);
    model.h5x          = zeros(nx,nx,nx,nx,nx,nx);
    model.hssxxx       = zeros(nx,nx,nx,nx);
    model.hsssxx       = zeros(nx,nx,nx);
    model.h4sx         = zeros(nx,nx);
    model.h5s          = zeros(nx,1);
end
model.eta              = eta;
model.vectorMom3       = auxOut.momEps;
model.ny               = ny;
model.nx               = nx;
model.mx               = auxOut.mx;
model.myx              = auxOut.myx;
model.ne               = ne;
model.params           = params;
model.transformY       = auxOut.transformY;    
model.transformX       = auxOut.transformX;    
model.Yss              = auxOut.Yss;
model.Xss              = auxOut.Xss;
model.orderApp         = orderApp;
if exist('selectY','var')
    model.derivs.gx        = derivs.gx(selectY,:);
    model.derivs.gxx       = derivs.gxx(selectY,:);
    model.derivs.gxxx      = derivs.gxxx(selectY,:);
    model.derivs.gxxxx     = derivs.gxxxx(selectY,:);
    model.derivs.gxxxxx    = derivs.gxxxxx(selectY,:);
    model.labely           = model.labely(selectY);
    model.g0               = model.g0(selectY,:);
    model.gx               = model.gx(selectY,:);
    model.gxx              = model.gxx(selectY,:,:);
    model.gxxx             = model.gxxx(selectY,:,:,:);
    model.g4x              = model.g4x(selectY,:,:,:,:);
    model.g5x              = model.g5x(selectY,:,:,:,:,:);
    model.gss              = model.gss(selectY,:);
    model.gssx             = model.gssx(selectY,:);
    model.gssxx            = model.gssxx(selectY,:,:);
    model.gssxxx           = model.gssxxx(selectY,:,:,:);
    model.gsss             = model.gsss(selectY,:);
    model.gsssx            = model.gsssx(selectY,:);
    model.gsssxx           = model.gsssxx(selectY,:,:);
    model.g4s              = model.g4s(selectY,:);
    model.g4sx             = model.g4sx(selectY,:);
    model.g5s              = model.g5s(selectY,:);
    model.transformY       = model.transformY(selectY);
    model.Yss              = model.Yss(selectY);
    model.ny               = length(model.Yss);
end
end


